%Main file to do empirical Kalman-filtering work in Leduc, Moran, and Vigfusson


ciCritValue = 6.25;  %90 percentile with 3 df. 


finalyear =2019;

%Inputs data files on current spot prices and futures prics.
%Spot prices can be extended here http://www.eia.gov/dnav/pet/pet_pri_spt_s1_d.htm
%Futures Prices.  12 and 24 Month Ahead Closing WTI Prices from NYMEX
%To extend, use Bloomberg CL12 and CL24.  



%% load in data. 
copyInData = 0;
 
% This section is for monthly spot prices. 
[sptdata,sptname,sptdatesstring]=tblread('spotwtileduc2020m4.prn',',');
sptdates = datenum(sptdatesstring,'mmm yy');
plot(sptdates,sptdata);

load wtibackdata.txt
wtibackdatadate = datenum(wtibackdata(:,1),wtibackdata(:,2),1);
plot(wtibackdatadate,wtibackdata(:,3),sptdates,sptdata)
datetick('x')
 
[wtidates,sptinx,histinx] = union(sptdates,wtibackdatadate);
plot(wtidates,[wtibackdata(histinx,3) ; sptdata(sptinx)])
wti = [wtibackdata(histinx,3) ; sptdata(sptinx)];

 
[aa,~,cc] = tblread('wtiforecastdatabest.csv','comma');
wtiForecast = aa;
forecastX= datenum(cc,'mm/dd/yyyy');


 %This section is for daily futures prices. 
[aa,bb,cc] = tblread('cloilfutleducf2020.csv','comma');
dcc = datenum(cc);
ixtwelve = contains(cellstr(bb),"CL12",'IgnoreCase',true);
oneyearfut = aa(:,ixtwelve);
oneyearfutDates = dcc;

ixtwentyfour = contains(cellstr(bb),"CL24",'IgnoreCase',true);
twoyearfut = aa(:,ixtwentyfour);
twoyearfutDates = dcc;


ixfront = endsWith(cellstr(bb),"CL1",'IgnoreCase',true);  %has to use endsWith to avoid CL11 and CL12
frontfut = aa(:,ixfront);

fut1dtz = dcc;
ftsprd = oneyearfut./frontfut;
[~ ,mftsprdnan,mfutsprddate] = bettermonthavg([frontfut  ftsprd],fut1dtz );

estdata = log(mftsprdnan);  
estdates = mfutsprddate;
ix = find(isnan(estdata(:,2)),1,'last');
estdata = estdata(ix+1:end,:);
estdates =estdates(ix+1:end);

%jmppoint = mod(3-month(estdates(1)),3); %to make sure that the first obs is end of quarter. 
% Matlab put the month function in the financial toolbox.  

yearmonthestdates = datevec(estdates(1));
jmppoint = mod(3-yearmonthestdates(2),3); %to make sure that the first obs is end of quarter. 


 %This will be the quarterly data used in making charts. 
 qestdata = [diff(estdata((jmppoint+1):3:end,1)) estdata((jmppoint+4):3:end,2)];
 qestdates =  estdates((jmppoint+4):3:end,:);

 

dix = datenum(sort(repmat(1989:finalyear,1,12))',repmat(1:12,1,length(1989:finalyear))',7);
inz = NaN + dix;
for zx= 1:length(dix)
    try
        inz(zx) = find( oneyearfutDates <dix(zx),1,'last');
        
        
    catch
    end
    
end

inz = inz(~isnan(inz));
ft1yr = log(oneyearfut(inz,:));

%some days will be missing a futures price.
%If so, then use the price from 1 to 5 days ago. 
for zx = 1:5
ftbackup = log(oneyearfut(inz-zx,:));
ft1yr(isnan(ft1yr)) =  ftbackup(isnan(ft1yr));
end
dix1yr = dix;


dix = datenum(sort(repmat(1989:finalyear,1,12))',repmat(1:12,1,length(1989:finalyear))',7);
inz = NaN + dix;
for zx= 1:length(dix)
    try
    inz(zx) = find( twoyearfutDates<dix(zx),1,'last');
    catch
    end
    
end

inz = inz(~isnan(inz));
ft2yr = log(twoyearfut(inz,:));


for zx = 1:5
ftbackup = log(twoyearfut(inz-zx,:));
ft2yr(isnan(ft2yr)) =  ftbackup(isnan(ft2yr));
end


dix2yr = dix;

[mprice,mpricenan,mdate] = bettermonthavg(wti,wtidates);

%yearstart  = find(month(mdate(1:15)) == 1,1,'first');
yearmonthmdate = datevec(mdate(1:15));
yearstart  = find(yearmonthmdate(:,2) == 1,1,'first');



qprice = mprice((yearstart+2):3:end,:);
qpricenan = mpricenan((yearstart+2):3:end,:);
qdate = mdate((yearstart+2):3:end);
 
qpriceavgnan =(1/3)*(mpricenan(yearstart:3:end-2)+mpricenan((yearstart+1):3:end-1) + mpricenan((yearstart+2):3:end));


 %Summary: We have constructed spot prices where we use the last month of each
 %quarter as our observed spot price.


yyavg    = log(qpriceavgnan(:,1))';
 
yy    = log(qpricenan(:,1))';
%Trim off missing observations from the front and end.
dq = qdate(find(~isnan(yy),1,'first'):find(~isnan(yy),1,'last'));

yy = yy(find(~isnan(yy),1,'first'):find(~isnan(yy),1,'last'));
if any(isnan(yy))
    error('missing prices');
end
 


% warning off for optimization
optz = optimset('fminunc');
optz = optimset(optz,'Display','off'); 
optz = optimset(optz,'Largescale','off'); 


%We use all oil prices since 1981 when Reagan deregulated the U.S. Oil
%Industry.  
zstart=  find(dq== datenum(1981,3,1));

stpoint = find(dq== datenum(1999,3,1));
enpoint = find(dq== datenum(2019,12,1));
 

sv1 = [0.0133    0.7126    0.0798    0.1556];
frbj = inline('[x(1:2) log(abs(x(3:4)))]','x');
revbj = inline('[x(1:2) exp(x(3:4))]','x');


%% Begin Estimation Part
% Code for estimation of univariate model
addpath code4rob

bestaltv = zeros(3,3,enpoint);

havebestanswer = 0;  %if you want to skip the estimation step. 
if havebestanswer == 1
    
    load havebestanswer 
else 
    %look for it.
    for zx = enpoint:-1:stpoint
        bjll  = @(x) -ll([0 x], diff(yy(zstart:zx)));
        sv = frbj(sv1);
        [bjbout,ll1,~,~,~,ftwo1]  = fminunc(bjll ,sv(2:end), optz);
        bestaltv(:,:,zx) =inv(ftwo1);
        bjbout = revbj([0 bjbout]);
        bjalt(zx,:) =  bjbout ;
    end
    
    
    
    
    
    
             
tsim = 10000;

%do you want to construct confidence intervals around the futures
%forecasts. If so set doconfidenceintervals = 1;  
 doconfidenceintervals = 1;  

 
 for zx = enpoint:-1:stpoint
     sv =frbj(bjalt(zx,:));
     
     %Here we filter the data to get the permanent and transitory
     %components.
     [llike, xi_tt_out] = ll(sv,  diff(yy(zstart:zx)));
     dp = diff(yy(zstart:zx));
     muplusvalt = dp(2:end)' -(xi_tt_out(2:end,1)-xi_tt_out(1:end-1,1));
     ix = find(dq(zstart+1:end) == datenum(1992,3,1));
     %we don't have a level for a starting point. So we set 1992Q3 equal to
     %log($20).
     percompalt = cumsum([log(20);muplusvalt(ix:end)]);
     bjmeanfut(zx,2) = percompalt(end)+(sv(2)^8)*xi_tt_out(end,1)+sv(1)*8;
     
     %Because there is a Jensen's inequality term we add in a simulation
     %step to get the expected value of exp(ft)
     bjmeanexpfut(zx,2) = mean(exp(bjmeanfut(zx,2)+sum( exp(sv(3))*randn(tsim,8)+ repmat((sv(2).^(0:7)),tsim,1).*exp(sv(4)).*randn(tsim,8),2)));
     
     bjmeanfut(zx,1) = percompalt(end)+(sv(2)^4)*xi_tt_out(end,1)+sv(1)*4;
     bjmeanexpfut(zx,1) = mean(exp(bjmeanfut(zx,1)+sum( exp(sv(3))*randn(tsim,4)+ repmat((sv(2).^(0:3)),tsim,1).*exp(sv(4)).*randn(tsim,4),2)));
     
     
     if doconfidenceintervals == 1
         
         %Now to do confidence intervals.
         %We pick parameter values from the normal distrbution for the parameter
         %values.
         origsv = frbj(bjalt(zx,:));
         bpicks = bsxfun(@plus,origsv(2:4),(chol(bestaltv(:,:,zx))'*randn(3,5000))');
         
         %only pick parameter values that are in the confidence interval.
         
         bpart = bsxfun(@minus,bpicks,origsv(2:4));
         part1 = sum((bpart*inv(squeeze(bestaltv(:,:,zx)))).*bpart,2);
         bcheck = bpicks(part1<ciCritValue,:);
         numtocheck = size(bcheck,1);
         maxsample = 500;
         if numtocheck > maxsample
             bcheck = bcheck(ceil(rand(maxsample,1)*numtocheck),:);
         end
         numtocheck = size(bcheck,1);
         
         
         altmeanexpfut2year = zeros(numtocheck,1);
         altmeanexpfut1year = zeros(numtocheck,1);
         
         ix = find(dq(zstart+1:end) == datenum(1992,3,1));
         dp = diff(yy(zstart:zx));
         parfor zp = 1:numtocheck
             
             
             sv =[0 bcheck(zp,:)];
             sv(2) = min(sv(2),0.99);
             sv(2) = max(sv(2),-0.99);
             [llike, xi_tt_out] = ll(sv,dp);
             muplusvalt = dp(2:end)' -(xi_tt_out(2:end,1)-xi_tt_out(1:end-1,1));
             percompalt = cumsum([log(20);muplusvalt(ix:end)]);
             altfutmean = percompalt(end)+(sv(2)^8)*xi_tt_out(end,1)+sv(1)*8;
             altmeanexpfut2year(zp) = mean(exp(altfutmean+sum( exp(sv(3))*randn(tsim,8)+ repmat((sv(2).^(0:7)),tsim,1).*exp(sv(4)).*randn(tsim,8),2)));
             altfutmean = percompalt(end)+(sv(2)^4)*xi_tt_out(end,1)+sv(1)*4;
             altmeanexpfut1year(zp) = mean(exp(altfutmean+sum( exp(sv(3))*randn(tsim,4)+ repmat((sv(2).^(0:3)),tsim,1).*exp(sv(4)).*randn(tsim,4),2)));
             
             
         end
         
         
         altmeanexpfut = [ altmeanexpfut1year altmeanexpfut2year];
         sa =sort(altmeanexpfut);
         
         maxaltgenfut(zx,1:2,1) =  sa(ceil(0.95*numtocheck),1:2);
         minaltgenfut(zx,1:2,1) = sa(ceil(0.05*numtocheck),1:2);
         maxaltgenfut(zx,1:2,2) =  sa(ceil(0.82*numtocheck),1:2);
         minaltgenfut(zx,1:2,2) = sa(ceil(0.17*numtocheck),1:2);
         meanaltgenfut(zx,1:2) = mean(sa(:,1:2));
     end
     
 end
 
 save havebestanswer
end


figure
textsize = 14;
mkfigtwo(dq(stpoint:enpoint),bjalt(stpoint:enpoint,:),bestaltv(:,:,stpoint:enpoint),textsize);

orient landscape
print -dpdf -bestfit finalfigures/Fig2

close all

yearticks = 2000:5:2020;
varp = bjalt(stpoint:enpoint,3).^2;
vart = bjalt(stpoint:enpoint,4).^2;
phi = bjalt(stpoint:enpoint,2);
vartempcomp = vart.*((1-bjalt(stpoint:enpoint,2)).^2)./(1-bjalt(stpoint:enpoint,2).^2);
volshare_base(stpoint:enpoint) = varp./(varp + 2*vart./(1+bjalt(stpoint:enpoint,2)));



figure
subplot(3,1,1)

plot(dq(stpoint:enpoint),volshare_base(stpoint:enpoint),'linewidth',3,'color',[94,204,243]/255)
set(gca,'fontsize',textsize)
axis([dq(stpoint) dq(enpoint)+365 0 1])  

set(gca,'xtick',datenum(yearticks,6,30), 'Layer', 'top')
datetick('x','keeplimits','keepticks')

bcov = bestaltv(:,:,stpoint:enpoint);
volshare_variance = zeros(length(varp),1);
for zip = 1:length(varp)
    vp = sqrt(varp(zip));  %exp(sqrt(varp(zip)));
    vt = sqrt(vart(zip)); %exp(sqrt(vart(zip)));
    zphi = phi(zip);
p1_dvolshare =  2*vp*vt/(vp+2*vt+zphi*vp); 
dvolshare = p1_dvolshare*[1 zphi+1  -(zphi+1)];

volshare_variance(zip) = dvolshare*squeeze(bcov(:,:,zip))*dvolshare';
end

hold on 
lwbnd = [dq(stpoint:enpoint) volshare_base(stpoint:enpoint)'-2*sqrt(volshare_variance)'];
upbnd = [dq(stpoint:enpoint) volshare_base(stpoint:enpoint)'+2*sqrt(volshare_variance)'];
cipatch = [lwbnd;upbnd(end:-1:1,:)];
patch(cipatch(:,1),cipatch(:,2),[0.85 0.85 0.85])
plot(dq(stpoint:enpoint),volshare_base(stpoint:enpoint),'linewidth',3,'color',[94,204,243]/255)
vaxis = axis;
axis([vaxis(1:2) -0.2 1])
plot(vaxis(1:2),[0 0],'k--')
orient landscape
orient landscape
print -dpdf -bestfit finalfigures/Fig3


figure
subplot(3,1,1:2)
fighandle = mkfigsix(dq(stpoint:enpoint),exp(yy(stpoint:enpoint)),dix2yr,exp(ft2yr),bjmeanexpfut(stpoint:enpoint,2),[minaltgenfut(stpoint:enpoint,2,1)  maxaltgenfut(stpoint:enpoint,2,1)]);

print -dpdf -bestfit finalfigures/Fig4v1

 figure
subplot(3,1,1:2)
fighandle = mkfigsixbw(dq(stpoint:enpoint),exp(yy(stpoint:enpoint)),dix2yr,exp(ft2yr),bjmeanexpfut(stpoint:enpoint,2),[minaltgenfut(stpoint:enpoint,2,1)  maxaltgenfut(stpoint:enpoint,2,1)]);
print -dpdf -bestfit finalfigures/Fig4v1

  
 RPstpoint = find(qestdates== datenum(1999,3,1));
 RPenpoint = find(qestdates== datenum(2019,12,1));

 clear svgood
svgood(RPstpoint:RPenpoint,1:5) =  [bjalt(stpoint:enpoint,:)  0*bjalt(stpoint:enpoint,1)];
conoptz = optimoptions('fmincon','Display','off'); 

 for zx = RPenpoint:-1:RPstpoint
     newsv = [0.8  0.7078    0.1049    0.1328 0.2];

    smpl = 1:zx;
    rvll =   @(params)llf_restat(qestdata(smpl,:)' ,params,[],[],4);
    [outvalues,llone(zx)] = fminunc(rvll,newsv , optz);
    svvalues(zx,:) =    outvalues;
 
    newsv = svgood(zx,:);
    smpl = 1:zx;
    rvll =   @(params)llf_restat(qestdata(smpl,:)' ,params,[],[],4);
    [outvalues,lltwo(zx)] = fminunc(rvll,newsv , optz);
    svvaluestwo(zx,:) =    outvalues;

    %LB = [-1 -1 0 0 0];
    %UB = [1 1 Inf Inf Inf];
    
    LB = [-1+eps -1+eps 0 0 0];
    UB = [1-eps 1-eps Inf Inf Inf];

    [outvalues,llfour(zx),~,~,~,~,storehess]  = fmincon(rvll,newsv ,[],[],[],[],LB,UB,[],conoptz);
    svvaluesfour(zx,:) =    outvalues;
    
     
     [tmpout,~,~,~,~,althess]  = fminunc(rvll ,outvalues, optz);
     samecoefbaseline(zx) = isequal(round(tmpout,4),round(outvalues,4));
     lastwarn('');
     
     
     strsderrbaseline(:,1,zx) =     sqrt(diag(inv(althess)));
    [warnMsg, ~] = lastwarn;
    if ~isempty(warnMsg)
        zx;
    ' bad constrained inverse';   
    end
     lastwarn('');
      
     strsderrbaseline(:,2,zx) =     sqrt(diag(inv(althess)));
      
         [warnMsg, warnId] = lastwarn;
   if ~isempty(warnMsg)
        zx;
    ' bad unconstrained inverse';  
    
      sqrt(diag(inv(althess(2:end-1,2:end-1))))
      zx;
      
    end
      
        
    
    
    
end
 
%These three points ended up at local maximums.  This gives them a chance
%to reconsider. 
for zx = [RPstpoint RPstpoint+2 RPstpoint+3  ]
    smpl = 1:zx;
    LB = [-1 -1 0 0 0];
    UB = [1 1 Inf Inf Inf];
    
    rvll =   @(params)llf_restat(qestdata(smpl,:)' ,params,[],[],4);
    [outvalues,lltest]  = fmincon(rvll,svvaluesfour(RPstpoint+1,:) ,[],[],[],[],LB,UB,[],conoptz);
    if lltest < llfour(zx)
        svvaluesfour(zx,:) =    outvalues;
        llfour(zx) = lltest;
        'fixed it';
        
        [tmpout,~,~,~,~,althess]  = fminunc(rvll ,outvalues, optz);
        strsderrbaseline(:,2,zx) =    sqrt(diag(inv(althess)))';
    end
end



%Constrained Models

 for zx = RPenpoint:-1:RPstpoint
    newsv = [0 0.5    0    0.1328 0.3];
    newsv = [newsv(1:2) newsv(4:5)];
    smpl = 1:zx;
    rvllnoperm =   @(params)llf_restat(qestdata(smpl,:)' ,[params(1:2) 0 params(3:4)],[],[],4);
  %  [outvalues,llnoperm(zx)] = fminunc(rvllnoperm,newsv , optz);
    LB = [-1 -1 0 0];
    UB = [1 1 Inf Inf];
     [outvalues,llnoperm(zx),~,~,~,~,storehess] = fmincon(rvllnoperm,newsv ,[],[],[],[],LB,UB,[],conoptz);
     outnoperm(zx,:) =    outvalues;
    
     [tmpout,llnopermunc(zx),~,~,~,althess]  = fminunc(rvllnoperm ,outvalues, optz);
     samecoefnoperm(zx) = isequal(round(tmpout,4),round(outvalues,4));
      strsderrnoperm(:,1:2,zx) =    [sqrt(diag(inv(storehess)))  sqrt(diag(inv(althess)))];

end



 for zx = RPenpoint:-1:RPstpoint
    newsv = [0 0.5    0.1    0.1328 0.3];
    newsv = newsv([1 3 5]);
    smpl = 1:zx;
    LB = [-1 0 0];
    UB = [1 Inf Inf];
    rvllnotemp =   @(params)llf_restat(qestdata(smpl,:)' ,[params(1) 0 params(2) 0 params(3)],[],[],4);
    [outvalues,llnotemp(zx),~,~,~,~,storehess] = fmincon(rvllnotemp,newsv ,[],[],[],[],LB,UB,[],conoptz);
    outnotemp(zx,:) =    outvalues;
     
     [tmpout,~,~,~,~,althess]  = fminunc(rvllnotemp ,outvalues, optz);
     samecoefnotemp(zx) = isequal(round(tmpout,4),round(outvalues,4));
     strsderrnotemp(:,1:2,zx) =    [sqrt(diag(inv(storehess)))  sqrt(diag(inv(althess)))];
    
    
end


  convalues(stpoint:enpoint,:)  = svvaluesfour( RPstpoint:RPenpoint,:);
   
  
sebaseline  =   squeeze(strsderrbaseline(:,2,end));
senotemp    =   squeeze(strsderrnotemp(:,2,end));
senoperm    =   squeeze(strsderrnoperm(:,2,end));

  
tabzout = [[0 0  outnotemp(end,[2 1 3])]' [outnoperm(end,[2 3]) 0  outnoperm(end,[1 4])]' convalues(end,[2 4 3 1 5])'];
%seout = [[0 0  senotemp([2 1 3])]' [senoperm([2 3]) 0  senoperm([1 4])]' sebaseline([2 4 3 1 5])];
seout = [[0 0  senotemp([2 1 3])']' [senoperm([2 3])' 0  senoperm([1 4])']' sebaseline([2 4 3 1 5])];

bigtabz = NaN(10,3);
bigtabz(1:2:10,:) = tabzout;
bigtabz(2:2:10,:)= seout;


rp_std = sqrt( tabzout(end,:).^2./(1-tabzout(end-1,:).^2));

sebaselinepass(stpoint:enpoint,:)  =   squeeze(strsderrbaseline(:,2,RPstpoint:RPenpoint))';


% SWP Code for Table. 
mkswtable([bigtabz;rp_std],'%0.4f')

 mkfigcoef(dq(stpoint:enpoint),convalues(stpoint:enpoint,:),sebaselinepass(stpoint:enpoint,:),textsize);

 orient landscape
print -dpdf -bestfit finalfigures/Fig7coef
 
  

close all
figure

yearticks = [2000 2010 2020];
subplot(2,2,1)
plot(dq(stpoint:enpoint),100*sqrt((1/1-convalues(stpoint:enpoint,1).^2).*(convalues(stpoint:enpoint,5).^2))./(convalues(stpoint:enpoint,3)),'k-','linewidth',1.25)
hold on
plot(dq(stpoint:enpoint),100*sqrt((1/1-convalues(stpoint:enpoint,1).^2).*(convalues(stpoint:enpoint,5).^2))./(convalues(stpoint:enpoint,4)),'r--','linewidth',1.75, 'color',[0.4 0.4 0.4]);
axis([dq(stpoint)  dq(enpoint)+365 0 6])
set(gca,'xtick',datenum(yearticks,6,30), 'Layer', 'top')
set(gca,'xticklabels',num2str(yearticks'))
datetick('x','keeplimits')
legend('Relative to $\sigma_{p}$','Relative to $\sigma_{\tau}$','location','Northeast','Interpreter','latex','fontsize',14)
ylabel('Percent','fontsize',14)
 %xlabel('Estimation End Point')
%title('\zeta')
 set(gca,'fontsize',14)
orient landscape
print -dpdf -bestfit finalfigures/FigRelVarBW

  
%figure
%mkfigcoeffourpanel(dq(stpoint:enpoint),convalues(stpoint:enpoint,:),sebaselinepass(stpoint:enpoint,:),textsize)%,bjalt(stpoint:enpoint,:));
%orient landscape
%print -dpdf -bestfit newfinalfigures/Fig8relvar



%% Survey Figure
  
figure

ix = sum(~isnan(wtiForecast),2);
wtiForecastNoZero = wtiForecast;
wtiForecastNoZero(isnan(wtiForecast)) = 0;
wtiForecast = sum(wtiForecastNoZero,2)./ix;
mwtiForecast= NaN+[wtiForecast wtiForecast];
mwtiForecast30= NaN+[wtiForecast wtiForecast];

for zq = 2:length(forecastX) 
     ix = find(and(forecastX<= forecastX(zq),forecastX>(forecastX(zq)-45)));
     if ~isempty(min(wtiForecast(ix)))

     mwtiForecast(zq,1) = min(wtiForecast(ix));
     mwtiForecast(zq,2) = max(wtiForecast(ix));
     end
    ix = find(and(forecastX<= forecastX(zq),forecastX>(forecastX(zq)-30)));
     if ~isempty(min(wtiForecast(ix)))

     mwtiForecast30(zq,1) = min(wtiForecast(ix));
     mwtiForecast30(zq,2) = max(wtiForecast(ix));
     end
     
 end


figure
subplot(3,1,1:2)
rangearea = [[forecastX; forecastX(end:-1:1)] [mwtiForecast30(:,1) ; mwtiForecast30(end:-1:1,2)]]  ;
 ix = ~any(isnan(rangearea'));
pp = patch(rangearea(ix,1),rangearea(ix,2),[0.9 0.9 0.9],'EdgeColor',[0.9 0.9 0.9]);
%plot(dix2yr,exp(ft1yr),'--','linewidth',2.5,'color',[241,65,36]/255)           % 1-yr ahead future
hold on
% I will provide this line
plot(dq(stpoint:enpoint)+30,bjmeanexpfut(stpoint:enpoint,1),':','color',[155 155 155]/255,'linewidth',2.5)   % 1-yr modeled
% end
hold on
vaxis = axis;
%plot(forecastX,wtiForecast,'o','color',[155 155 155]/255,'markerfacecolor',[155 155 155]/255,'markersize',4);
%axis([forecastX(10)  dq(enpoint)+60 vaxis(3) 180])  
datetick('x','keeplimits');
legend('Survey Range', '1-year Ahead Futures Prices (Model)','Location','southwest')

%legend('Survey Range','1-year Ahead Futures Prices (Model)', 'Location','northeast')
orient landscape
ylabel('Dollars Per Barrel','Fontsize',16)
set(gca,'Fontsize',16)
shortyearticks = 2006:2:2020;
set(gca,'xtick',datenum(shortyearticks,6,30))
axis([forecastX(1) datenum(2020,12,30) 0 160])
datetick('x','keeplimits','keepticks');
box on
orient landscape
print -dpdf -bestfit finalfigures/newSurveybw 
  
  
  
  
  